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The epidemic-type aftershock sequence model (ETAS) is a simple stochastic process modeling seismicity, 
based on the two best-established empirical laws, the Omori law (power law decay ~ l/f^^" of seismicity 
after an earthquake) and Gutenberg-Richter law (power law distribution of earthquake energies). In order to 
describe also the space distribution of seismicity, we use in addition a power law distribution ~ l/r^^*^ of 
distances between triggered and triggering earthquakes. The ETAS model has been studied for the last two 
decades to model real seismicity catalogs and to obtain short-term probabilistic forecasts. Here, we present an 
exact mapping between the ETAS model and a class of CTRW (continuous time random walk) models, based 
on the identification of their corresponding Master equations. This mapping allows us to use the wealth of 
results previously obtained on anomalous diffusion of CTRW. After translating into the relevant variable for 
the ETAS model, we provide a classification of the different regimes of diffusion of seismic activity triggered 
by a mainshock. Specifically, we derive the relation between the average distance between aftershocks and the 
mainshock as a function of the time from the mainshock and of the joint probability distribution of the times 
and locations of the aftershocks. The different regimes are fully characterized by the two exponents 6 and /i. 
Our predictions are checked by careful numerical simulations. We stress the distinction between the "bare" 
Omori law describing the seismic rate activated directly by a mainshock and the "renormalized" Omori law 
taking into account all possible cascades from mainshocks to aftershocks of aftershock of aftershock, and so on. 
In particular, we predict that seismic diffusion or sub-diffusion occurs and should be observable only when the 
observed Omori exponent is less than 1, because this signals the operation of the renormalization of the bare 
Omori law, also at the origin of seismic diffusion in the ETAS model. We present new predictions and insights 
provided by the ETAS to CTRW mapping that suggest novel ways for studying seismic catalogs. Finally, we 
discuss the present evidence for our predicted sub-diffusion of seismicity triggered by a main shock, stressing 
the caveats and limitations of previous empirical works. 



I. INTRODUCTION 



The spatio-temporal complexity of earthquakes is often invoked as an illustration of the phenomenon of critical 
self-organization with scale-invariant properties JUSSBEIl- This concept points to the importance of developing 
a system approach in which large scale properties can emerge from the repeating interactions occurring at smaller 
scales. Such ideas are implemented in models proposing links between the physics of earthquakes and concepts 
of statistical physics, such as critical points, self-organized criticality, spinodal decomposition, critical depinning, 
etc., in order to explain the most solidly established facts in the phenomenology of earthquakes, of which we cite 
the three most important. 

• LAW 1 : The Gutenberg-Richer law |0] states that the cumulative distribution of earthquake magnitudes 
m sampled over broad regions and large time intervals is proportional to lO"*"", with a fo-value h k, \. 
Translating into energies E with the correspondence m — (2/3) logig E-\- constant leads to a power law 
- l/E^ with B« 2/3. 

• LAW 2: Omori law for aftershocks fT] states that the rate of earthquakes triggered by a mainshock decays 
with time according to an inverse power l/i^ of time with an exponent p sa 1. 

• LAW 3: The earthquakes are clustered in space along hierarchical fault structures |8] and their spatial 
distribution over long times can be approximately described by a fractal dimension close to 2.2 (in three 
dimensions) ||9|]. 
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There are many other empmcal "laws" but these three characterize the very fundamentals of seismicity in size, 
time and space. 

We should immediately point out that these three laws come with significant caveats. 

1 . There have been on-going controversies on the universality of the exponent B or 6-value of the Gutenberg- 
Richterlaw [10.. 11.,]. 

2. The exponent p of Omori law exhibits a large variability from one aftershock sequence to another aftershock 
sequence and is found typically in the range from 0.3 to 2. We note however that not all these values, 
especially the extreme ones, automatically reflect a bona-fide power law decay and one should exert caution 
in attributing too much confidence to them. 

3. The view that geological faults and earthquake hypocenters are fractal objects is now recognized to be a 
naive description of a much more complex reality in which a hierarchy of scales occur with possibly different 
organizations at different scales |8]. 

In addition, a major difficulty for making progress in modeling and predicting earthquakes is that these three 
and other laws may be "explained" by a large variety of models, with many different mechanisms. For instance, 
with respect to the first two laws, we observe the following. 

• There are many mechanisms that create a power law distribution of earthquake sizes (see for instance the list 
of mechanisms described in chapter 14 of Ref. jl^ ). 

• Omori law is essentially a slowly decaying "propagator" describing a long-time memory of past events 
impacting on the future seismic activity. Such slow power law time decay of the Omori propagator may 
result from several and not necessarily exclusive mechanisms (see iflBil and references therein): pore-pressure 
changes due to pore-fluid flows coupled with stress variations, slow redistribution of stress by aseismic creep, 
rate-and-state dependent friction within faults, coupling between the viscoelastic lower crust and the brittle 
upper crust, stress-assisted micro-crack corrosion il4ll5lil . slow tectonic driving of a hierarchical geometry 
with avalanche relaxation dynamics 1 16|, etc. 

The zeroth order description of earthquakes is to consider a single isolated homogeneous fault on which earth- 
quakes are recurrent to accommodate the long-term slow tectonic loading. But faults are not isolated and the most 
conspicuous observation is that earthquakes interact and influence each other on complex fault structures. Un- 
derstanding these interactions is essential for understanding earthquakes and fault self-organization. However, the 
full impact of interactions between earthquakes is still far from being well understood. The simplest and clearest 
observation of earthquake interaction is provided by aftershocks whose phenomenology is captured by Omori law 
(LAW 2). Indeed, aftershocks are the most obvious and striking signature of the clustering of the seismicity in 
time and space, and are observed after all large shallow earthquakes. Most aftershocks are triggered a few hours 
or days after the mainshock. However, due to the very slow power law decay of the rate of aftershock, known as 
the Omori law |7], aftershocks can be triggered up to a hundred years after the mainshock 1 17]. Aftershocks often 
occur near the rupture zone of the mainshock with a variety of focal mechanisms suggesting that they are actually 
on separate structures They are also sometimes triggered at very large distances from the mainshock 

j20, 21, 22, 23, 24]. As an example. Hill et al. ]20] observed aftershocks of the Landers earthquake as far as 
1250 km from the epicenter Similarly to the temporal distribution of aftershocks, a power-law distribution seems 
to describe well the distribution of distances between pairs of events 122] . Since a power-law decays slowly, it 
describes a slow decay of the probability of observing aftershocks at large distances to the mainshock. 

Thus, Omori law can be considered as the simplest and best established description of earthquake interactions 
of a certain kind. The question we investigate is whether it can be used fruitfully to explain a larger variety of 
earthquake interactions beyond the class of observations that were used to establish it. In a series of papers 
l^,*??], we find that Omori law for aftershocks plus the constrain that aftershocks are distributed according to the 
Gutenberg-Richter power law for earthquake size distribution independently of the magnitude of their progenitor 
is enough to derive many of the other empirical "laws," as well the variability of the p-exponent. Here, we test the 
potential of this approach to account for and to quantify observations on aftershock diffusion. 

Aftershock diffusion refers to the phenomenon of expansion or migration of aftershock zone with time 12^12^ 
|30, 31, 32, 33, 34, 35, 36]. Immediately after the mainshock occurrence, most aftershocks are located close to the 
rupture plane of the mainshock, then aftershocks seem to migrate away from the mainshock, at velocity ranging 
from 1 km/hour to 1 km/vea r i36ll37il . Note that this expansion is not universally observed, but is more important 
in some areas than in others 13111 . 
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The diffusion of aftershocks is usually interpreted as a diffusion of the stress induced by the mainshock, either 
by a viscous relaxation process |37], or due to fluid transfer in the crust |35, 38, 39]. Another interpretation of 
the expansion of aftershocks is given by Dieterich |40J, who reproduces the Omori law decay of aftershocks and 
the expansion of the aftershock zone with time, using a rate and state friction law and assuming that the rate of 
aftershocks is proportional to the stress rate. In his model, the expansion of aftershock zone arises from the non- 
uniform stress induced by the mainshock. Another alternative explanation is that the diffusion of aftershocks is 
mainly due to the occurrence of large aftershocks, and to the localization of secondary aftershock close to the 
largest aftershocks, as observed by Ouchi |34|. The apparent diffusion of the seismicity may thus result from a 
cascade process; the mainshock triggers aftershocks that in turn trigger their own aftershocks, and thus lead to an 
expansion of the aftershock zone. 

In the present paper, we investigate the epidemic time aftershock sequence model (ETAS), and show that the 
cascade of secondary aftershocks can indeed explain the reported diffusion of aftershocks. The ETAS model was 
introduced by Kagan and Knopoff 141] (in a slightly different form than used here) and Ogata 1 42 1 to describe the 
temporal and spatial clustering of seismicity. This model provides a tool for understanding the clustering of the 
seismic activity, without arbitrary distinction between aftershocks, foreshocks and mainshocks. In this model, all 
earthquakes are assumed to be simultaneously mainshocks, aftershocks and possibly foreshocks. Each earthquake 
generates aftershocks that decay with time according to Omori law, which will in turn generate their own after- 
shocks. The seismicity rate at any given time and location is given by the superposition of aftershock sequences of 
all events impacting that region at that time according to space-time "propagators." The additional ingredient in the 
version of the ETAS model that we study is that the number of aftershocks per earthquake increases exponentially 
cx 10"™ with the magnitude m of the mainshock (i.e., as a power law oc E'^"/^ of the energy released by the 
mainshock), in agreement with the observations |43','44'1. Since the energy of an earthquake is a power law of its 
rupture length, this law expresses the very reasonable idea that the number of events related to a given earthquake 
is proportional to a power of its volume of influence. The value of the exponent a controls the nature of the seismic 
activity, that is, the relative role of small compared to large earthquakes. Few studies have measured a in seismicity 
data ll43ll45ll46ll . This parameter a is often found close to b |43] or fixed arbitrarily equal to b |41, 47]. In the case 
where a is close to the Gutenberg-Richter 6-value, this law also reproduces |47] the self-similar empirical Bath' s 
law 1 48], which states that the average difference mjv/ — niA in size between a mainshock and its largest aftershock 
is 1.2 magnitude units, regardless of the mainshock magnitude: niA — tum — 1-2. If a < 6, small earthquakes, 
taken together, trigger more aftershocks than larger earthquakes. In contrast, large earthquakes dominate earth- 
quake triggering \f a>b. This case a > b has been studied analytically in the framework of the ETAS model by 
Ref. 1I27I1 and has been shown to eventually lead to a finite time singularity of the seismicity rate. This explosive 
regime cannot however describe a stationary seismic activity. 

A natural way to tame this singular behavior is to introduce an upper cut-off for the magnitude distribution 
at large magnitudes, mirroring the cut-off toq used for the low-magnitude range. The physical argument for 
introducing this cut-off is based on the finiteness of the maximum earthquake that the earth is capable of carrying. 
The specific way of introducing such a cut-off (abrupt or smooth with a transition to a power law with larger 
exponent or to an exponential taper) is not very important qualitatively because all these laws will regularize 
the singular behavior and make the average branching ratio finite. Such regularization with a maximum upper 
magnitude then allows a > b. The special case a = b required for Bath's law to hold exactly can not therefore be 
excluded. 

However, based on a recent re-analysis of seismic catalogs using the powerful collapse technique, one of us 14^ 
has presented strong evidence that a is strictly smaller than 6. In this paper, we will therefore consider only the 
case a < b and take a = 0.5 specifically in our numerical simulations. In this regime a < b, Bath' s law cannot be 
reproduced because the average difference in size between a mainshock and its largest aftershock increases with 
the mainshock magnitude. For a < fe, it is easy to show that Bath's law is replaced by tjia = {a/b)mM— constant, 
where toj\/ and niA are the magnitudes of the mainshock and of the largest aftershock. Tests of this prediction will 
be reported in a future publication but we expect that distinguishing this modified Bath's law from Bath's law will 
be a difficult task due to the limited range of the studied magnitudes as well as the dependence of the distribution 
of niM ~ ™a on the magnitude thresholds chosen for the mainshocks and for the aftershocks [49jl . 

We assume that the distribution of all earthquakes follow the Gutenberg-Richter distribution and take this dis- 
tribution of aftershock sizes to be independent of the magnitude of the mainshock. Therefore, an earthquake can 
trigger a larger earthquake, albeit with a small probability. This model can thus describe a priori both aftershock 
and foreshock sequences. The ETAS model has been calibrated to real seismicity catalogs to retrieve its parameters 

]Ei m E3 113 El m m |54] and to give short-term probabilistic forecast of seismic activity by extrapolating 

past seismicity into the future via the use of its space-time propagator L4lil55ll56ll . 
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The ETAS model is a branching model which exhibits different regimes EfiT depending upon the value of the 
branching ratio n, defined by the average number of primary aftershocks per earthquake. The critical case n = 1 
corresponds to exactly one primary aftershock per earthquake, when averaging over all mainshock magnitudes 
larger than a threshold mo. Let us stress that n is an average quantity which does not reflect adequately the 
large variability of the number of aftershocks per main shock, as a function of its magnitude. Indeed, the number 
of aftershocks per mainshock increases exponentially fast as a function of the mainshock magnitude, so that large 
mainshocks will have significantly more than n aftershocks. For a = 0.5, a magnitude 7-earthquake gives typically 
10 times more direct aftershocks than a magnitude 5, and 100 times more direct aftershocks than a magnitude 3- 
earthquake. The increase in triggered seismic activity with the magnitude of the mainshock is obviously stronger 
for a larger value of a. Note that these numbers refer to aftershocks of the first generation; the total number of 
triggered events is larger by the factor 1/(1 — n) ~ 10 (for n w 0.9 which is typical), due to the cascades of 
secondary aftershocks. Notwithstanding this large variability, the average number n of primary aftershocks per 
earthquake controls the global regime. For n exactly equal to 1, seismicity is at the border between death and 
growth. In the sub-critical regime n < 1, since each earthquake triggers on average less that one aftershock, 
starting from a large event, the seismicity will decrease with time and finally die out. The super-critical n > 1 
corresponds to more that one primary aftershock per earthquake on average. Starting from a large earthquake, after 
a transient regime, the average seismicity will finally increase exponentially with time [26,1 . but there is still a finite 
probability for aftershock sequences to die out. 

The numerical simulations reported below have been performed with a = 0.5. It is probable that a good fit 
to seismic data is obtained by using a value of a sa 0.8 larger that the value 0.5, as reviewed and documented 
recently by one of us 146]. We have checked that results similar to those presented below hold true qualitatively for 
larger values 0.5 < a < 1. Such larger values of a lead however to stronger fluctuations which are more difficult 
to handle numerically because the variance of the number p{m) of direct triggered aftershocks defined below in 
(|3} becomes undefined for a > 0.5. A full understanding of this regime requires a special treatment that will be 
reported elsewhere. 

Sornette and Sornette j25|l studied analytically a particular case of this model, without magnitude and spatial 
dependence, and they considered only the subcritical regime n < 1. Starting with one event at time t = 
and considering that each earthquake generates an aftershock sequence with a "local" Omori exponent p = 1 + 
9, where 6 > 0, they studied the decay law of the "global" aftershock sequence, composed of all secondary 
aftershock sequences, i.e., by taking into account that the primary aftershocks can create secondary aftershocks 
which themselves may trigger tertiary aftershocks and so on. They found that the global aftershock rate decays 
according to an Omori law with an exponent p =1 — 6'<1, uptoa characteristic time ll25lE<3l 



nr(i 



(1) 



|l-n| 

and then recovers the local Omori exponent p = 1 + 6 for time larger than t*. Helmstetter and Sornette |2^ 
extended their analysis to the general ETAS model with magnitude dependence, and considered both the sub- 
and the super-critical regime, but still restricted the analysis to the temporal distribution of the seismicity, without 
spatial dependence. In the sub-critical regime, they recovered the crossover found by Sornette and Sornette jzsll . In 
addition, Helmstetter and Sornette fl^ give the explicit mathematical formula for the gradual transition between 
the Omori law with exponent p = 1 — 9 for t ^ t* to the Omori law with exponent p = 1 + 9 for t ^ t*. 
This smooth transition can be observed in figure |2] on the line calculated for t* = 10^ days with n < 1. t* can 
thus be viewed as the time where the apparent exponent p of the Omori law is approximately in between the two 
asymptotic values 1 — 9 and 1 + 6'. A more rigorous mathematical definition 1 26] is that t* is the characteristic time 
scale such that /3t* is the dimensionless variable of the Laplace transform (with variable f3) of the seismicity rate. 

In the super-critical regime, Helmstetter and Sornette |26] found a novel transition between a power-law decay 
with exponent p = 1 — 6' at early times, similar to the sub-critical regime, to an exponential increase of the 
seismicity at large times. The regime where a > 6 or equivalently 2a/3 > B has been found to lead to a new 
kind of critical stochastic finite-time-singularity iZTll . relying on the interplay between long-memory and extreme 
fluctuations. Recall that the number of aftershocks per earthquake increases as a power law oc iJ^"/-^ of the energy 
released by the mainshock whereas the number of earthquakes of energy E decreases as the Gutenberg-Richter law 
cx . Intuitively, when 2a/3 > B, the increase in the rate of creation of aftershocks with the mainshock 

energy more than compensate the decrease of the probability to get a large mainshock when the mainshock energy 
increases. This theory based solely on the ETAS model has been found to account for the main observations (power 
law acceleration and discrete scale invariant structure) of critical rupture of heterogeneous materials, of the largest 
sequence of starquakes ever attributed to a neutron star as well as of some earthquake sequences IIStIi . 
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In the sequel, we extend the analytical study of the temporal ETAS model to the spatio-temporal 

domain. To model the spatial distribution of aftershocks, we assume that the distance between a mainshock and 
each of its direct aftershock is drawn from a given distribution, independently of the magnitude of the mainshock 
and of the delay between the mainshock and its aftershocks. For illustration but without loss of generality for 
the mapping to the continuous time random walks model (CTRW) discussed later, we shall take a power law 
distribution of distances between earthquakes. We take the simplest and most parsimonious hypothesis that space, 
time and magnitude are decoupled in the earthquake propagator Our first result is to establish a correspondence 
between the ETAS model and the CTRW, first introduced by Montroll and Weiss f57\ and used to model many 
physical processes. We then build on this analogy to derive the joint probability distribution of the times and 
locations of aftershocks. We show analytically that, for sufficiently short times t < t*, the average distance 
between a mainshock and its aftershock increases subdiffusively as i? , where the exponent H depends on 
the local Omori exponent 1 + 9 and on the distribution of the distances between an earthquake and its aftershocks. 
We also demonstrate that the local Omori law is not universal, but varies as a function of the distance from the 
mainshock. Due to the diffusion of aftershocks with time, the decay of aftershock is faster close to the mainshock 
than at large distances. These non-trivial space-time couplings occur notwithstanding the decoupling between 
space, time and magnitude in the "bare" propagator, and is due to the existence of cascades of aftershocks. 

A recent work of Krishnamurthy et al. L58.1 substantiates the general modeling strategy used here of representing 
the space-time dynamics of earthquakes by an effective stochastic process (the ETAS model) entirely defined by 
two exponents (corresponding to our /i and H{6,ii) defined below), where fi is the exponent of the power law 
distribution of jumps between successive active sites and H is the (sub-)diffusion exponent. Indeed, Krishnamurthy 
et al. 1 58] show that the Bak and Sneppen model and the Sneppen model of extremal dynamics (corresponding 
to a certain class of self-organized critical behavior fl2!|) can be completely characterized by a suitable stochastic 
process called "Linear fractional stable motion." Beyond recovering the scaling exponents of this model, the 
stochastic process strategy predicts the conditional probabilities of successive activations at different sites and thus 
offers novel insights. We note that this approach with the Linear fractional stable motion is extremely close in 
spirit as well as in form to our approach mapping the ETAS model to the CTRW model. The ETAS model can thus 
be taken to represent an effective stochastic process of the complex self-organization of seismicity. 

II. THE ETAS MODEL 

A. Definitions and specific parameterization of the ETAS model 

We assume that a given event (the "mother") of magnitude rrii occurring at time ti and position gives birth to 
other events ("daughters") of any possible magnitude chosen with some independent Gutenberg-Richter distribu- 
tion at a later time between t and t + dt and at point rztdr to within dr at the rate 

(t>m,{t-U,f-ri) = p{mi) -^{t^ti) $(f-r*) . (2) 

We will refer to (prm {t — ti, r — Ti) both as the seismic rate induced by a single mother or as the "bare propagator". 
It is the product of three independent contributions: 

1. p{mi) gives the number of daughters born from a mother with magnitude rrii. This term will in general be 
chosen to account for the fact that large earthquakes have many more triggered events than small earthquakes. 
Specifically, we take 

p{mi) = K io«(™^-'"») , (3) 

which, as we said earlier, is justified by the power law dependence of the volume of stress perturbation as a 
function of the earthquake size, a quantifies how fast the average number of daughters per mother increases 
with the magnitude of the mother. 

2. ^'(i — ti) is a normalized waiting time distribution giving the rate of daughters born at time t ~ ti after the 
mother. The normalization condition reads ^ dt ^[t) = \. 5'(< — ti)dt can thus be interpreted as the 
probability for a daughter to be born between t and t + dt from the mother that was born at time ti. '^{t — ti) 
embodies Omori law; it is the "bare" or "direct" Omori law 



(4) 
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where 9 > Q and H{t) is the Heaviside function. 

3. $(r — fi) is a normahzed spatial "jump" distribution from the mother to each of her daughter, quantifying 
the probabiUty for a daughter to be triggered at a distance |r — | from the mother. Specifically, we take 

$(f) = ^ , (5) 

which has the form of an (isotropic) elastic Green function dependence describing the stress transfer in an 
elastic upper crust. The exponent /i is left adjustable to account for heterogeneity and the possible complex 
modes of stress transfers. The normalization condition reads / dr $(r) = 1 where the integral is carried out 
over the whole space. 

The physical justification for this decoupled model (|3 in which 0m, (i ^ ti,r — fi) is the product of three 
independent distributions is that elastic waves propagate at kilometers per second and thus almost instantaneously 
reset the stress field after a large main shock. In other words, there is a well-defined separation of time scales 
between the time of propagation of seismic waves (seconds to minutes) which control the convergence to a new 
mechanical equilibrium after the main shock and the time scales involved in aftershock sequences (hours, days, 
months to many years). The spatial dependence in ^ reflects the stress redistribution. This new stress field 
then relaxes slowly and more or less independently from point to point leading to the local Omori law '^{t ~ ti). 
Notwithstanding this argument, the decoupling in (|3 between the local responses in magnitudes, space and time 
is mostly performed because of its simplicity. It constitutes an approximation that should be checked and relaxed 
in future studies. 

We assume a distribution P{m) of earthquake sizes expressed in magnitudes m which follows the Gutenberg- 
Richter distribution 

P{m) = b In(lO) I0-b{m-mo) ^ (g) 

with a 6- value usuaUy close to 1. mo is a lower bound magnitude below which no daughter is triggered. 



B. The branching ratio n 

A key parameter of the ETAS is the average number n of daughter-earthquakes created per mother-event, 
summed over all possible magnitudes. As we shall see, it is also natural to call it the "branching ratio". To 
see this, consider the integral of the seismic rate {t — ti,r — fi) induced by one earthquake over all times after 
ti, over all spatial positions and over all magnitudes > ttiq, which must give by definition the average number 
n of direct (or primary) daughter-earthquakes created per mother-event independently of its magnitude. For a < b 
and using (|2j, (|3} and (|6}, it is exactly given by 

df dt I drrii P{mi) 4>mi{t - ti,f- f) = / drrii P{mi) p{mi) = , (7) 

Jt, J mo J mo O-a 

since the two integrals over time and space contribute each a factor 1 by the normalization of ^ and This result 
is identical to that found in absence of spatial dependence of (f>mi{t — ti) with respect to f — f due to the 
factorization of the rate p, time and space $ dependences l26ll . The branching ratio has also been evaluated in 
the case where the magnitude distribution follows a gamma distribution fs?]. 

We stress again that n is an average quantity which does not reflect the large fluctuations in the number of 
aftershocks from events to events. Indeed, large events with magnitudes M produce in general many more after- 
shocks than small events with magnitude m < M, simply because p{M) ^ p{iti) if Af > m (see the exponential 
dependence (jSjl of p{m) on the magnitude to). 



C. Numerical simulation of the spatial ETAS model 



The ETAS model has been simulated numerically using the algorithm described in Refs. js^E^]. Starting with 
a large event of magnitude M at time t = 0, events are then simulated sequentially. At any given time t, we 
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calculate the conditional seismic rate X{t) defined by 

where K = n{b — a)/b, and ti and nii are the times and magnitudes of all preceding events that occurred at time 
ti < t. Note that we use the bare propagator because the sum in is performed exhaustively on the complete 
catalog of past events. The time of the following event is then determined according to the non-stationary Poisson 
process of conditional intensity \{t), and its magnitude is chosen according to the Gutenberg-Richter distribution 
with parameter h. To determine the position in space of this new event, we first choose its mother randomly among 
all preceding events with a probability proportional to their rate of aftershocks (pm . {t — ti) evaluated at the time of 
the new event. Once the mother has been chosen, we generate the distance r between the new earthquake and its 
mother according to the power-law distribution "I>(r) given by (|5|i. The location of the new event is determined by 
assuming an isotropic distribution of aftershocks. By this rule, it is clear that new events tend to be close in general 
to the last large earthquakes, leading to space clustering. 

Note that this two-steps procedure is equivalent to but more convenient for a numerical implementation than the 
one-step method, consisting of calculating at each point on a fine space-covering grid the seismic rate, equal to the 
sum over all preceding mothers weighted by the bare space ^(r) and time ^'(t) propagators given by Q and (|3J; 
after normalizing, these rates then provide to each grid point a probability for the event to occur on that point. The 
equivalence between our two-step procedure and the direct calculation of the seismic rates is based on the law of 
conditional probabilities: probability of next event (A) = probability of next event conditioned on its mother (event 
B) X probabiHty of choosing the mother, i.e., P{A, B) = P{A\B) x P{B). 

Figure[2shows the result of a numerical simulation of the ETAS model which exhibits a diffusion of the seismic 
activity. We simulate a sequence of aftershocks and secondary aftershocks starting from a mainshock of magnitude 
A/ — 7, with the following parameters: 9 = 0.2, 6 = 1, a = 0.5, n — \ and /i = 1. At early times, aftershocks are 
localized close to the mainshock, and then diffuse and cluster close to the largest aftershocks. This (sub-)diffusion 
is extremely slow, as we shall quantify in the sequel. Our purpose is to provide a theory for this process based on 
the ETAS model. This theory will be tested by numerical simulations. 

The different regimes are illustrated by Figure |2] which shows the seismicity rate N{t) for the temporal ETAS 
model studied by |25, 26] obtained by summing the seismic activity over all space for the 3 cases n < 1 (sub- 
critical), n — 1 (critical) and n > 1 (super-critical). The sub-critical regime is characterized by the existence of 
the time scale t* given by Q. There is no difference between the critical case n = 1 and the sub-critical case for 
t < t* (see figure|3. Indeed, the difference between the sub-critical regime and the critical regime can be observed 
only for t > t*. A simple way to see this is to realize that the critical regime n = 1 gives t* ~ +oo, meaning that, 
in the critical regime, one is always in the situation t < t*. 

It is interesting to note that the spatial distribution of epicenters shown in the right panel of figure ^ has the 
visual appearance of a fractal set of points. This is confirmed by the calculation of the correlation dimension of 
this set of = 3000 points generated in the time interval [30, 70] yrs, which is found approximately equal to 
D2 = 1.5 ± 0.05 over more than two decades in spatial scales, as shown in figure|3] If we use instead all 30, 000 
events of the simulation performed up to time t = 70 yrs, we find D2 — 1 .85±0.05 while the correlation dimension 
of the geometrical set made of the epicenters of the 10, 000 last events (time interval [7, 70] yrs) is D ^ 1.7± 0.05, 
also over more than two decades in scale. These values are similar to those reported for 2D maps of active fault 
systems i8[[6l| l62ll63ll . and are in good agreement with D2 values in the range [1.65, 1.95] measured for aftershocks 
epicenters [641. The fractal clustering of the earthquake epicenters, according to the ETAS model, occurs because 
of a self-similar process taking place on many different scales. However, the description of this multi-scale process 
solely in terms of a single fractal dimension fails to fully embody the complex spatial superposition of local 
"singularities" associated with each aftershock on the one hand and finite-size effects (stemming from the finite 
lifetime of each aftershock sequence) on the other hand. Each event indeed creates its cloud of direct aftershocks 
which can be characterized by its singular exponent 1 — /i for /i < 1 and for /i > 1, defined by the scaling 
oc rdr/r^^'^ oc of the "mass" of the cloud with its radius R. Finite-size effects and randomness have 

been documented to generate realistic but sometimes spurious fractal signatures [65.. 66.. ■67„.68J . This problem 
requires a special study which is left for another work. 



8 



D. Relationship witli tiie space-independent ETAS model 

The spatial ETAS model reduces to the space-independent ETAS model solved in 1 26] by integrating the dressed 
propagator obtained below over all space. In the Fourier representation (see expression illi ). this corresponds to 
putting the wavenumber k to zero. Indeed, for fc = 0, the Fourier transform amounts to perform a simple integration 
over all space. Since <l>(fc = 0) = 1, expression i2H derived below reduces to the form studied at length in \2§]. 
Therefore, all results reported previously hold also for the version of the space-dependent ETAS model studied 
here, when averaging over the whole space. This is an important property that all the solutions discussed below 
must obey. 

in. MAPPING OF THE ETAS MODEL ON THE CTRW MODEL 

In order to study the space-time properties of the ETAS model, it is very useful to use an exact correspondence 
between the ETAS model and the continuous time random walk (CTRW) that we establish here. In this way, 
we can adapt and use the wealth of results previously derived for the CTRW. But first, let us demonstrate the 
correspondence between the ETAS and CTRW models. For this, our strategy is to derive the Master equation for 
both models and show that they are identical. 

A. The Master equation of the ETAS model 

The ETAS model can be rephrased by defining the rate _,,„(t — ti,r — fi) at which a given event (the 
"mother") of magnitude rrii > niQ occurring at time ti and position fi gives birth to other events ("daughters") of 
specified magnitude m at a later time between t and t + dt and at point f to within an infinitesimal volume |dr|. 
Note that the only difference with respect to the previous definition Q is that we now specify also the magnitude 
m of the daughter. (j>mi-tm{t ~ ti,f — fi) is given by 

(t)m,^rn{t - ti,f~ f) = p{m, m) *(t - ti) $(f - fi) , (9) 

where ^(i — ti) and $(f — ri) are the same as previously while 

p{mi m) — P{m) p{mi) . (10) 
With the parameterization Q and (|5Jl, this reads 

p{m, ^ m) = n In(lO) (6 - a) 10°("'-"") iQ-Hm-mo) 

Let us consider the case where there is an origin of time i = at which we start recording the rate of earthquakes, 
assuming that a large earthquake has just occurred at i = and somehow reset the clock. In the following 
calculation, we will forget about the effect of events at times prior to t = and count all aftershocks that are 
created only by this main shock. 

Let us call Nm{t, f)dt dm dr the number of earthquakes occurring between t and t + dt of magnitude between m 
and 171 + dm inside of box of volume \df\ centered at point f. Nm{t, f) is the solution of a self-consistency equation 
that formalizes mathematically the following process : an earthquake may trigger aftershocks; these aftershocks 
may trigger their own aftershocks, and so on. The rate of seismicity at a given time t and position r is the result of 
this cascade process. The self-consistency equation that sums up this cascade reads 

Nmit,f)=S{t,f,m)+ f dr f dm' f dr (l>rn'^m{t - T,f - f) N^,{T,f) . (12) 
J J niQ J 

The rate Nm{t, f) at time t and position r is the sum over all induced rates from all earthquakes of all possible 
magnitudes that occurred at all previous times and locations propagated to the present time t and to the position f 
of observation by the corresponding bare propagator. The induced rate of events per earthquake that occurred at an 
earlier time t and position f is equal to 4>m'^m{t — t, f— f). The source term S{t, f) is the main shock plus the 
background seismicity, if any. In absence of background seismicity, a main earthquake which occurs at the origin 
of time < = at position r = with magnitude M gives 



S{t, r, m) = 6{t) 5{m - M) 5{f) 



(13) 
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where 5 is the Dirac distribution. Other arbitrary source functions can be chosen. 

The source term corresponding to a single mainshock is indeed the deha function M3\ rather than the direct 
Omori law created by this mainshock in direct lineage. To see this, notice that the direct Omori law is recovered 
from M2\ by replacing Nm'{T,f') in the integral by S{t,r,m) given by il3\ . This shows that the difference 
between the renormalized and the direct Omori laws comes from taking into account the secondary, tertiary, etc., 
cascade of aftershocks. 

As we have seen, a key assumption of the ETAS model is that the daughters born from a given mother have 
their magnitude drawn independently of the magnitude of the mother and of the process that give them birth, 
with a probability given by the Gutenberg-Richter distribution (|6jl. The consequences resulting from relaxing this 
hypothesis will be reported elsewhere. Keeping this assumption, it can be shown i69ll that for a < &/2 an ensemble 
of realizations will obey 

Nmit, r) = P(m) N{t, f) , for i > , (14) 

which makes explicit the separation of the magnitude from the time and space variables. N{t, r) is the number 
of events at position r at time t of any possible magnitude. Expression MA\ means that the Gutenberg-Richter 
distribution is preserved at all times. That J14I I holds for the ETAS model stems from the fact that the waiting time 
^'(i) distribution and jump size $(r) distribution (|5} are independent of the magnitudes and that fluctuations 
in the seismicity rate are not too wild for a > 6/2. Note that, in a more complex model in which time, space 
and magnitudes are interdependent, expression il4\ would become a mean-field approximation, in which the 
fluctuations of the rates induced by the fluctuations of the realized magnitudes of the daughters factorize from 
the process. 

Putting il4l in il2\ . for t > when the source term S{t, r, m) is identically zero, one can simplify by P{m) 
and obtain 

N{t,r) = J dr' dr (f>{t-T,f-r') N{T,r') , t>0, (15) 

where 

/•oo 

(p{t - T,r- f') ^ / dm' P{m')(j)^'{t-T,r-r') . (16) 

J mo 

Equation il5\ is nothing but the expectation (or statistical average, i.e., average over an ensemble of realizations) 
of expression (|8jl, with the definition N{t, r) = E[A(t) Therefore, the Master equation obtained here gives 

us only the first moment of the space-time dynamics of seismicity. It is not difficult to derive the equations for the 
variance and covariance of the seismic rate as well as higher moments. 

The value of the source term alt — Q that should be incorporated in M5\ requires more care. Indeed, a naive 
treatment would give a source term d{t)S{m — M)5{r) / P{M) obtained by simply dividing by P{m), expressed 
at TO i\f due to the Dirac distribution 5{m — M). However, this source term still depends on m via the Dirac 
distribution 5{ra — M) and is thus unsuitable as a source term of the equation ( I15> which is independent of to. 
In order to circumvent this difficulty, one has to get rid of the Dirac distribution 5{m — M). The corresponding 
procedure has been described in details in Ref. i23l and consists in applying the integral operator dm 4>{j3, r) 

to J12I I. where 4){(3, r) is the Laplace transform with respect to the time variable of (f){t, r). In this way, the Dirac 
distribution 5{m — M) is regularized. Identifying with the results of Ref. L261 . we obtain that N{t, f) is solution 
of il5\ with a source term 

SM{t,r) = 6ir)S{t)p{M)/n, (17) 

where p{M) is defined in (|3jl and n is given by Q. Thus, the complete Master equation for the number N{t, f) of 
events at position f at time t of any possible magnitude is solution of 

N{t,r}=^SM(t,r) + Jdf' J^dT(f>{t-T,r-f')N{T,r'), t>0, (18) 

N{t, f) is the "dressed" or "renormalized" propagator, obtained by summing the bare Omori propagator over all 
possible aftershock cascades. N{t, f) can also be called the renormalized Omori law 125]. 

The essential assumption used to derive ( I12t is that the fluctuations of the earthquake magnitudes in a given 
sequence can be considered to be decoupled from those of the seismic rate. This approximation can be shown to 
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be valid for a < 6/2 i69ll . for which the random variable p{mi) has a finite variance. In this case, any coupling 
between the fluctuations of the earthquake energies and the instantaneous seismic rate provides only sub-dominant 
corrections to the equation (I12> . For a > 6/2, the variance of p{nii) is mathematically infinite or undefined as 
p{mi) is distributed according to a power law with exponent h/a < 2. In this case, the Master equation (I12> is 
not completely correct as an additional term must be included to account for the effect of the dependence between 
the fluctuations of earthquake magnitudes and the instantaneous seismic rate. Our results are presented below for 
a = 0.5 which belongs to the first regime a < 6/2. For a > 6/2, Ref. has shown that the renormalization of 
the bare propagator into the dressed propagator is weaker than for a < 6/2, all the more so as a — > 6. Preliminary 
numerical simulations for a > b/2 shows that our results presented below hold qualitatively but with a reduction 
of the observed spatial diffusion exponent compared to the value predicted from the Master equation approach 
developed here. This regime a > 6/2 is probably relevant to the real seismicity i43l l45l l46ll . even if a precise 
estimation of a is very difficult. 



B. A Master equation of the CTRW model 

We now demonstrate that the self-consistent mean field equation (I18> is identical to the Master equation of a 
continuous-time random walk (CTRW). Random walks underlie many physical processes and are often the basis 
of first-order description of natural processes. The CTRW model, which is a generalization of the naive model of 
a random walker which jumps by ±1 spatial step on a discrete lattice at each time step, was introduced by 
and investigated by many other workers The CTRW considers a continuous distribution of 

spatial steps as well as time steps (which can be seen either as waiting times between steps or as durations of the 
steps). The CTRW model is thus based on the idea that the length of a given jump, as well as the waiting time 
Ti = ti — ti^i elapsing between two successive jumps are drawn from a joint probability density function (pdf) 
t), which is usually referred to as the jump pdf. From a mathematical point of view, a CTRW is a process 
subordinated to random walks under the operational time defined by the process {ti}. 

From 0(r, t), the jump length pdf ^{r) = dt <j){r, t) and the waiting time pdf ^(t) = f df (j){r, t) can 
be deduced. Thus, $(r)(ir produces the probability for a jump length in the interval {f,r + df) and 'if{t)dt the 
probability for a waiting time in the interval {t, t + dt). When the jump length and waiting time are independent 
random variables, this corresponds to the decoupled form (l>{r,t) — If both are coupled, a jump of 

a certain length involves a time cost or, vice versa in a given time span the walker can only travel a maximum 
distance. With these definitions, a CTRW process can be described through a Master equation (see I74i l75ll76il for 
a review and references therein) which turns out to be given by an equation which is identical to jl8> . 

This connection between the ETAS model of earthquakes and a model of random walks provides an important 
advance for the understanding of spatio-temporal earthquake processes, as it allows one to borrow for the deep 
knowledge accumulated in past decades on random walks. In the same spirit, polymer physics acquired its status 
as a fundamental physical problem from its previous status of an applied field of research in chemistry when Flory, 
Edwards, de Gennes, des Cloizeaux and others showed how to formulate problems in polymer physics in the 
language of random walks and how to extract novel results. In the sequel of this article, we use this analogy to 
provide a wealth of new predictions as well as new questions for earthquake aftershocks. 

In the context of the CTRW, we have the following correspondence. 

• N{t, r) is the pdf for the random walker to just arrive at position r at time t. 

• The source term Si\i{t, r) given by (I17> denotes the initial condition of the random walk, here chosen to be 
at the origin of space at time t — Q. The constant p{M)/n adds the possibility via the parameter M to have 
more than one initial walker at the origin. 

• In the CTRW context, the Master equation ( I18> states that the pdf N{t, r) of just having arrived at position r 
at time t comes from all possible paths in number N{T,f') having crossed a position at an earlier time r, 
weighted by a transfer or propagator function (f){t — r, r — r*) describing all the possible steps of the random 
walker from (r, f*) to {t, f). 

It is important to stress that N{t, r) defined above is different from the standard quantity W{t, r) usually studied 
in random walk problems, defined as the probability to find the random walk at position f'at time t. The relationship 
between N{t, r) and W{t, r) is 





pt~t' 




f dt' 


1 - / dt" 'fit") 


N{t',r) 


JO 
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The term 1 — /q * dt" '^{t") in bracket is the probability for the walker not to jump in the time interval [t' , t] 
and the integral in the right-hand-side of ( I19> means that the probability W{t, r) for the random walker to be at 
position r at time t is the sum over all possible scenarios in which the walker just arrives at r at an earlier time t' 
and then does not jump until time t. In the context of earthquake aftershocks, W{t, r) is the probability that an 
event at r has occurred at a time t' <t and that the whole system has remained quiescent from t' to t. 
In the Fourier-Laplace domain (see below), expression M9\ reads 

W{P,k)^ -^N{p,k). (20) 

In general, the CTRW models transport phenomena in any heterogeneous media. It has for instance been used 
successfully for describing the behavior of chemical species as they migrate through porous media 
insight, it is rather natural that it can be applied to the "transport of stress" through the heterogeneous crust and 
thus to the description of the anomalous diffusion of seismic activity. 

Table|I]synthesizes the correspondence between the ETAS and CTRW models and then draws its consequences. 



C. Experimental verifications of the cross-over between the two power law Omori decays in photoconductivity in 
amorphous semi-conductors and in fractal stream chemistry using the correspondence between the ETAS and CTRW 

model 

The crossover from an Omori law for t < t* io for t > t* found in fzsl, with t* g iven 

by Q has actually a counterpart in the CTRW. This behavior was first studied by Scher and Montroll ITlll in a 
CTRW with absorbing boundary condition to model photoconductivity in amorphous semi-conductors As2Se3 and 
an organic compound TNF-PVK finding 6' w 0.5 and 9 = 0.8 respectively. In a semiconductor experiment, electric 
holes are injected near a positive electrode and then transported to a negative electrode where they are absorbed. 
The transient current follows exactly the transition for t < t* to for t > t* found for Omori law 

for earthquake aftershocks in the ETAS model. In the semiconductor context, the finiteness of t* results from the 
existence of a force applied to the holes while in the ETAS model it results from a finite distance 1 — n to the 
critical point n = 1 in the subcritical regime. When the force goes to zero or n ^ 1, t* — > +oo. 

A similar transition has been recently proposed to model long-term time series measurements of chloride, a 
natural passive tracer, in rainfall and runoff in catchments |80]. The quantity analogous to the dressed Omori 
propagator is the effective travel time distribution h{t) which governs the global lag time between injection of the 
tracer through rainfall and outflow to the stream. h(t) has been shown to have a power-law form h{t) ~ 
with m between -0.3 and 0.2 for different time series |81]. This variability may be due to the transition between 
an exponent 1 — 6* at short times to 1 + 6* at long times 1,80,1 . where 9 is the exponent of the bare distribution of 
individual transition times. 



D. General and formal solution of the spatial ETAS model 

Let us solve ( I18l l for the number N{t, r) of events at position r at time t of any possible magnitude. Recall that 
N{t, r) can also be interpreted as the dressed Omori propagator. Extending J26il to the spatial domain and also in 
analogy with the standard approach to solve the CTRW, the Laplace-in-time Fourier-in-space transform iV(/3, k) 
of N{t, r) is given by 

N{prk)-^4^^, (21) 

1 - 71*(/3)$(fc) 

where SM{f3,k) is the Laplace Fourier transform of the source SM{t,r) given by ilH and '^{(3) (respectively 
^(fc)) is the Laplace (respectively Fourier) transforms of (respectively "I>(r)). For a mainshock of magnitude 
M occurring at time t = and position r — 0, the source term is thus Sm{P, k) = p{M)/n. The only difference 
between expression M\\ and the Laplace-Fourier transform of the pdf of the CTRW of just having arrived at r at 
time t occurs when the branching ratio n is different from 1. In general, solutions of CTRW models are expressed 
for n — 1 and for the variable W{t^ r) which is simply related to N{t, r) according to M9\ . Using (I19> and ( I21> 
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leads to 

/? 1 - 71*(/3)$(fc) 

In the following, we exploit (I22t to obtain analytical solutions of the spatial ETAS model in different regimes, 
that provide specific predictions on the conditions necessary for observing aftershock diffusion. In addition, we 
provide specific predictions on the exponent H of the diffusion law R ^ that are tested by numerical simula- 
tions. 



IV. CRITICAL REGIME n = 1 
A. Classification of tiie different regimes 

Numerous works on the CTRW have investigated many possible forms for and 3>(r) and have provided the 
asymptotic long time and large scale dependence of W{t, r) (see j 74ll7^ l76l ItsIi and references therein). Here, 
we restrict our discussion to the cases where both '^{t) and $(r) have power law tails as given by Q and (|5}- The 
long-time and large scale behavior of the ETAS and CTRW are controlled by the behavior of the Laplace-Fourier 
transforms for small (3 and small 

Two cases must be distinguished depending on the exponent /i controlling the weight of the tail of ^{r). 

• For /i > 2, the variance ((r)^) — of the jump size distribution exists. To leading order in fc = |A;|, "I>(A;) 
can be expanded as 

^{k)^l-(j'^k^ +0{k°) , witho>2. (23) 

• For /i < 2, the variance ((r)^) is infinite. This regime of "long jumps" leads to so-called Levy flights. In 
this case, to leading order in fc = |fc|, $(fc) can be expanded as 

<l(fc) = l-(T^fc^ + 0(fc°) , where0</i<2, with o > /i , (24) 
where cr is a characteristic distance defined by 

d[T{l-^i)Yl>^, o<M<i, 

r(A<-l) sin(7r/x/2) ' -L < M < ^ ■ 

For a distribution '^{t) of waiting times of the form of a local Omori law with exponent <1, \t'(/3) can be 
expanded for small /3 as 

= 1 - {(3c'Y + Oip"^) , with uj>l. (26) 

1/9 

where cJ is proportional to c up to a numerical constant c' — c (r(l — 0)) in the case 6 < 1. 

Putting the leading terms of the expansions of $(fc) for small |fc| and of ^{(3) for small (3 in (I21t gives 

1 — n + nypc)" + na^k^^ 
The corresponding VK(/3, k) is obtained from J22t by 

W{(3, k) = SmW, k) ^7.,.. ^.^ . (28) 

The critical regime n = 1 gets rid of the constant term 1 — n in the denominator of (|^} and ( I28> . This case is 
analyzed in details below. 
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The regime n ^ \ introduces a characteristic time t* given by Q. In the sub-critical regime, equation (I27> can 
be rewritten as 



where r* is defined by 



" {\^n) 1 + {I3t*)0 + {kr*Y ■ ^^^^ 



For t < t* and ?■ < r*, the dressed propagator is given by the same expression as for the critical case and all our 
results below hold. For large times t > t* and large distances r > r*, we can factorize ( I29l l as a product of a 
function of time and a function of space 

^~ {1 + {f3t*)0) {1 + {kr*y) ■ ^ ' 

Thus, there is no diffusion in the sub-critical regime for t > t* and r > r*. We shall not analyze further this trivial 
regime n < 1 and t > t* and will only analyze the case i < t* . If there is the need, the cross-over can be calculated 
explicitly using (|27|i. 

In order to get the leading behavior of N(t, r) from that of W{t^ r), we see from (I21> and M2\ that N{j3^ k) — 
WiP, k) « /Ji-^c'-^ W{(3, k). The inverse Laplace transform of l//3^ is l/[r(6l) Using the fact 

that the Laplace transform of df /dt is (3 times the Laplace transform of f{t) minus /(O), we get N{t^ r) as the 
derivative of a convolution 

In ( I32t . we have dropped the Dirac function coming from the inverse Laplace transform of the constant term /(O), 
which provides a contribution only at the origin of time t — Q. Note that the operator ^ dt' ^^^,y-e is 
nothing but the so-called fractional Riemann-Liouville derivative operator of order \ — 9 applied to the function 
W{t, r) of time t and is usually denoted QD\^^W{t, r). 



B. The standard diffusion case 9 > 1 and ^ > 2 

The standard diffusion process is recovered for 9 >1 (for which the average waiting time is finite) and for /i > 2 
(for which the variance of the jump length is finite). In this case, iV(/3, k) = -0^^^- For an impulsive source 
leading to Sm{P, k) = constant, this is the Laplace-Fourier transform of the standard diffusion propagator 

« cxp[- (f)Vi^i] , where D = a^c' , (33) 

where d is here the space dimension. This solution is valid for \ir\/\/Di not too large. For larger values, large 
deviations lead to corrections with the power law tail of the input jump distribution $(r) ~ defined 
in along the lines presented for instance in fl2il (section 3.5). This regime is not relevant to the aftershock 
problem for which usually Q < 9 < 1. 

C. Long waiting times {9 < 1) and finite variance of the jump sizes > 2) 
Putting the leading terms of the expansions of $(fc) i23\ and of ^(/3) i26i in illi gives 

^(^'^^) = WTW '''' 
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The expression (I34> can be inverted with respect to the Fourier transform, and then inverted with respect to the 
Laplace transform using Fox functions SEl. The solution for W{t, f) in one dimension is given for instance in 
L76il in terms of an infinite sum 



k 



k=0 



e{k + i)/2) 



where 



(35) 



(36) 



andD = CT/c'«/2. 

Expression ( I35t and many others below involve the Gamma function of negative arguments. We recall that 
the Gamma function T{u) can be analytically continued to the whole complex plane, except for the simple poles 
u — 0, —1, —2, —3, ... Thus, r(u) is defined everywhere but at these poles. In order to get the expression of 
the Gamma function for negative arguments, one can use two formulae: r(l — u) x r(u) = tt/ sin(7ru) and 
r(l + u) = ur{u). Both these formulae are valid for all points with the possible exception of the arguments at 
poles 0,-1, -2,... For instance, r{-9) = T{l-e)/{~_^ = -[7r/6lsin(7r6')]/r(6'), for < 6 < 1. 

Expression ( I35> can be rewritten as a Fox-function 18311 



W(t,z) 



1 1 rrl,0 

2D if 



(1-0/2,0/2) 
(0,1) 



(37) 



whose asymptotic dependence for large z, obtained from a standard theorem of the Fox function (equation (1.6.3) 
of ^), 



W{t,z) 




(38) 



is in agreement with the result of Roman and Alemany f84l and Barkai et al. |82] for a space dimension dj — 1, 
including the dependence in the power law prefactor to the exponential. The exponential dependence W{t, r) ~ 

exp 1^— const {r / Dt^^^)^^^ in ( I38> holds in arbitrary dimensions df, the only modification occurring in the 

prefactor whose power of z change with the space dimension d/ as JsT 



WdAt,z) 



exp 



(39) 



The expression of N{t, r) can be obtained from W{t, f) using the fractional Riemann-Liouville derivation J32> 
of order 1 — 9. Inserting expression i35i in ( I32> and using the expression of the fractional Riemann-Liouville 
derivative operator o_D" applied to an arbitrary power t^, i.e., oD^f^ 



r(i+A') 
r(i+p-Q) 



" , we obtain 



N{t,r) 



(-1) 



2Dt^-l ^^klT{{l-k)0/2) 



(40) 



Expression ( I40> can be used to evaluate N{t, r) for small z, but the numerical evaluation of ( I40l l is impossible for 
larg e z. In order to obtain the asymptotic behavior of N{t^ r), expression ( I40> can be rewritten as a Fox-function 



Nit,r) 



2Dt^ 



H 
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{9/2,6/2) 
(0,1) 



Employing again the standard theorem of the Fox function (equation (1.6.3) of 
N{t, r) for large distances r such that r > Dt^/'^ is given by 



(41) 

), the asymptotic behavior of 
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The exponential dependence N{t, r) ^ exp ^— const {r/Dt^/^) 2-e j in (I42> holds in arbitrary dimensions. 

This expression becomes incorrect for very large distances because it would predict an exponential or slightly 
super-exponential decay with r. This cannot be true as the global law cannot decay faster than the local law (|5}- 
The reason for (142 > to become incorrect at large distances is that the expansion of N{(3, k) for small |fc| (large 
distances) given by \2>A\ has been truncated at the order k^. There is however a subdominant term cx fc'^ that 
describes the power law tail of the local law (|5|l and also of the global law asymptotically. A similar situation 
occurs in the application of the central limit theorem for sums of N random variables with po wer law d istributions 
with exponents /i > 2 fT2il : the distribution of the sum 5 is a Gaussian in its bulk for \S\ < \/N\nN and crosses 
over to a power law with tail exponent ^ for larger S. In a similar way, the cross-over of N{t, r) to the asymptotic 
local power law ^ can be recovered by an analysis including the subleading correction oc k'^ to the expansion 

(El- 
Expression ( I40> shows that the global rate of seismicity cannot be factorized as a product of a distribution of 
times and a distribution of distances. This space-time coupling implies that the seismic activity diffuses with time, 
and that the decay of the rate of aftershocks depends on the distance from the first mainshock. This coupling of 
space and time stems from the cascade of aftershocks, from the primary aftershocks to the secondary aftershocks 
to the tertiary aftershocks and so on. 

FigureHpresents the decay of the seismic activity N{r, t) obtained using expression ( I40l i for small z and expres- 
sion (l42t for large z, as a function of the time from the mainshock and as a function of the distances r. Close to the 
mainshock epicenter, expression ( I40I I predicts that the global seismicity rate decays with time as the renormalized 
Omori law 

(43) 

The same decay is found at any fixed point r for times t > {\f\/D)^^^. At all times, the same decay l/t^^^^^ 
is also obtained by measuring the aftershock seismicity in a local box at a distance from the main shock origin 
increasing with time as r ^ (this is nothing but putting z = constant in ( I40» . At large distances r > Dt^^^, 
the global decay law is different from a power-law decay. Figure |4] shows that the rate of aftershocks presents 
a truncation at early times, which increases as the distance r increases. At large times, the rate of aftershocks 
recovers the power-law decay (I43> . We stress that a fit of the global law N{r, t) over the whole time 

interval by an Omori law would yield an apparent exponent p < \ — 6/2 that decreases with r. 
Integrating (|^J over the whole one-dimensional space, we recover the global Omori law 

N{t) = j drN{t,r) ^ ^ (44) 

found in lEslEfill . Thus, we have found an additional source of variability of the exponent p of the Omori law: if 
measured over the whole catalog, we should measure p = 1 — 6* in the critical regime n ~ 1 while p = 1 — 9/2 
is slightly larger when measured in certain time- and space-windows, as described above. Thus, in this regime, 
pruning of catalogs may lead to continuous change from the value 1 — 9 to 1 — 9/2. In addition, as we have 
mentioned, the cross-over in time may lead to still smaller apparent exponents, thus enhancing the impression of 
variability of the exponent p. In reality, this range of p-values are seen to result from the complex spatio-temporal 
organization of the aftershock seismicity of the ETAS model. These results should lead us to be cautious when 
analyzing real catalogs with respect to the conditions and regimes under which the analysis is performed. 

There is another observable that characterizes how an aftershock sequence invades space as a function of time. 
Expression (I40> indeed predicts a sub-diffusion process quantified by 

(|rf ) ^ t^^ , (45) 

with H = 9/2 since the natural variable is z given by (I36> . Indeed, expression ( I4Q> tells us that, up to a global 
rescaling function of time, the rate of aftershocks is identical for a fixed value of z. Thus, any aftershock structure 
diffuses according to ( 145 1 . 

This prediction is checked in Figure|5]by numerical simulations. 1000 synthetic catalogs have been generated 
with fj, = i, 9 — 0.2 and n — 1. The average distance between the first mainshock and its aftershocks as a 
function of the time from the mainshock has been averaged over these 1000 simulations. The theoretical diffusion 
exponent is H = 9 /2 = 0.1, in good agreement with the asymptotic behavior observed in the numerical simulation. 
In practice, in order to minimize the effect of fluctuations and optimize the speed of convergence, we estimate 
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numerically exp[(ln which is also expected to scale as exp[(ln ^ t^/"^ due to the simple scaling form of 

ED. 

This problem has also been solved exactly in I85ll in the context of the so-called fractional Fokker-Planck 
equation, which amounts to replace the distribution ^{r) of jumps (|5|i by a Gaussian function. This fractional 
Fokker-Planck equation allows one to introduce the possibility of bias or drift in the CTRW and therefore in the 
aftershock sequence. 



D. Exponential waiting time distribution and long jump size Levy distribution (/i < 2) 

This case with an exponential distribution 

^il{t)^\ e~^* (46) 

of waiting times with a Levy distribution <I>(r) = _Lp(|r|) of jump sizes with tail exponent /i < 2 has been 
investigated by Budde et al. 1861 . One finds 

(|r=l2)i/2 _ ^^/^^ ^ (47) 

corresponding to a superdiffusion regime with Hurst exponent H = > 1/2. The full distribution function 
W{t, r) corresponding to the critical regime n = 1 is known for At >> 1: 



(48) 



The corresponding N{t, r) is obtained from ( I20> . The Laplace transform of the exponential distribution ( I46t is 

= A/(/3 + A). We thus get 

N{(3,k)^{(3 + \) W{p,k) , (49) 

and thus 

^ ^ dW{t r) ^ ^ ^ _ ^^^^ 
ot 

Expression J50> together with ( I48> predicts a diffusion law r ^ with H = which is in good agreement 
with our simulations. At large times |r| ^ (At)^/^, N{t,r) w A W{t,'r) ^ given an apparent local 

Omori exponent 9 = 1 — 1 //i. This offers a new mechanism for generating Omori law for aftershocks from purely 
exponential local relaxation but with a heavy distribution of jump sizes. This power-law decay should be observed 
only at a fixed distance r or over a limited domain from the mainshock in the regime of large times. 

Integrating over the whole space, J dr W{t, r) = 1 which gives N{t) = S{t) + A equal to a constant seismic 
rate. This results from an initial mainshock at t = leading to the cascade of aftershocks adjusting delicately to 
this constant rate for the critical value n = 1 of the branching parameter In the sub-critical regime n < 1, the 
Omori law integrated over space gives instead N{t) cx exp[— (1 — n)Xt], showing that the characteristic decay time 
1/(1 — n) A of the dressed Omori propagator N{t) becomes much larger (much longer memory) that the decay 
time 1/ A of the bare Omori propagator. 

For /i > 2, we recover the standard diffusion corresponding to 9 > 1 and /i > 2 discussed in section lTV Bl 



E. Long waiting times (9 < 1) and long jump sizes (Levy flight regime for /i < 2) 



Putting the leading terms of the expansions of $(fc) and of ^(/3) in Hit gives 



Niprk) = SMk) (^,,). ■ (51) 



The corresponding W{/3, k) is given by 
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Equation i52\ has been studied extensively in the context of the CTRW model as a long wavelength | fc | — > and 
long time /3 approximation to investigate the long time behavior of the CTRW. Kotulski [87] has developed 
a rigorous approach, based on limit theorems, to classify the asymptotic behaviors of different type of CTRWs 
and justifies the approximation i52\ for the long time behavior. Barkai |88] has studied the quality of the long 
wavelength |fc| ^ and long time /3 — > approximation i52\ by solving the exact CTRW problem for the 
case when the waiting time distribution ^'(t) is a one-sided stable Levy law of index 9 with the same tail as 
(0} and the distribution of jumps is a symmetric stable Levy of index /i with the same tail as (|5}. Their 
Laplace and Fourier transforms, that appear in the denominator of i22\ . are respectively ^(/3) = exp[— /3^] and 
= exp[— |fc|^/2]. Note that the long wavelength |fc| — > and long time (3^0 approximation gives 
1 — exp[— (c'/3)^] exp[— |(Tfc|^] — (c'/3)^ + |(Tfc|^, which recovers i5H . By comparing the exact solution of 
(I21> for ^'(t) and $(f^ of the above Levy form with that of the long wavelength |fc| — > and long time /? — > 
approximation i52\ . Barkai fsS*] finds that certain solutions of (I52> diverge on the origin, a behavior not found 
for the corresponding solutions of (12 U . In addition, certain solutions of the full equation (12 U converge only very 
slowly for /i < 1 to the solutions of the long-time approximation (I52> . These results validate our use of the 
asymptotic long time behavior with respect to the scaling laws but provide a note of caution if one needs more 
precise non-asymptotic information. In this case, such information can be obtained by a suitable analysis of the 
full equation i2l\ . 

Using power counting, expression i52\ predicts a diffusion process ( 145 > with exponent 

9 

H ^ - . (53) 

M 

This prediction is checked by numerical simulation of the ETAS model in the critical regime n = 1, with 9 — 0.2, 
1-1 = 0.9, shown in figure|6l The average distance between the first mainshock and its aftershocks as a function of 
the time from the mainshock indeed increases according to i45i with an exponent H in very good agreement with 
the prediction H ~ 9 / ji — 0.22. As the form of the denominator in \52\ is independent of the space dimension, 
the prediction ( I53l l is vaUd in any space dimension. 

The natural variable for the expansions given below allowing to compute N{t, r) is 

z = , (54) 

n 



where D = a/c"^/f' and c' = c (r(l - 9)) 



i/e 



1. z-expansion of the solution 
W{t, r) can be obtained as the following sum (equation (5.10) of jsojll 



H-oo 



Wit,f) = 4^ V (-1)™ S^^^ COS 



ni—0 

Applying i32\ to ( I55> term by term in the sum, we get 



T{m9 + 1) 



-(m^ + 1) 



N{t,r) 



J2i-irz'- 



m—O 



r(m^ + i) 

T{{m + 1)9) 



cos 



-(TO/i+ 1) 



The asymptotics 



T{mfi + fi + 1) T{m9 + 1) T{mfi + + 1) T{{m + 1)9) 



r{m9 + 9 + l)r(m^+l) 



T{{m + 2)9)r(mn+l) 



(55) 



(56) 



(57) 



show that the series i55\ and ( I56l l exist only for fi < 9. It can be shown that these series exist for all z in this 
case. This series converges very slowly for large z but the Fade summation method f9^ can be used to improve 
the convergence of (I56> in the case (i < 9, and can also be used to estimate (I56> in the case fj, > 9 for which the 
series diverges. 
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The space integral J dr N{t, r) over the whole one-dimensional volume V, with N{t^ r) given by (in, recovers 
the global Omori law 

J^drN{t,r)^ . (58) 

Note the non-trivial phenomenon in which the superposition of all aftershock activities transforms the local Omori 
law or "bare propagator" Q ^'(t) ^ ^ttf i'^'^o the global Omori law or "dressed propagator" -pk^- This effects 
was predicted in 1 25, 26] in the version of the ETAS model without space dependence. These results are consistent 
with the claim of section 111 Dl according to which all results reported previously for the version of the ETAS 
model without space dependence hold also for the version of the space-dependent ETAS model studied here, when 
averaging over the whole space. 

The asymptotic behavior for {r] ^ D t~ (i.e., z ^ 1) and jj. < 9 is obtained by keeping only the first non-zero 
term (m — 1) in i56\ which is convergent for all z in the case fj. < 6 

At fixed large | f] and for i < \ f/D\'^ , this predicts a local Omori law with exponent p = 1 — 29. 



2. 1/ z-expansion of the solution 



We use the theory of Fox functions f83T to obtain N{t, r) as an infinite series in 1/ z. For this, we first rewrite 
expression ( I56> as a Fox function L83.1 



N{t,r) 



R H. 



'2,2 



t/2 



(1^,1/m),(1,1) 



(60) 



where R{z) indicates the real part of z. 

The 1/ z expansion of N{t, r) can be obtained using the dual expansion of the Fox function ( I60> (expression 
(3.7.2) of JH) 



m— 



fl Z 



r(l - (to + 1)^) sin((TO + l)/X7r/2) 
r(-TO6l) 



TT cos(TO7r/2) 



(61) 



m! sin((TO -h l)7r/At) V(Q - (m + 1)61/^) 

This expansion exists only for [i> Q (conditions of page 71 below eq. (3.7.2) of ll83ll '). This is easily checked by 
the behavior of an asymptotics similar to ( I57> . Note that the series (16 U is not defined in the special case /i = 1 
due to the presence of the ill-defined ratio r(0)/r(0) and a different approach is required, such as the integral 
representation of VF(t, r) developed in |89]. The global Omori law obtained by integrating over the whole space 
( 16 1> is again N{t) ~ \ jt^~^ as expected from the analysis of the ETAS model without space dependence |26]. 
Keeping only the largest term of ( 16 H for large z, we obtain the asymptotic behavior for small distances r < 

D 



N{t,r) 



r(l - 2/i) sin(7r/i) sin(7r6') r(l 



e'er TT^ 



1 



da ^iT{6 - sin(7r/^) (i/c')i~''+^/^ 



for ^ < 0.5 
for 0.5 < /i < 2 



(62) 



Note that for r < D i^/^ and 0.5 < /i < 2, the leading behavior of N{t, r) is independent of r. 
Equation \62\ thus predicts an apparent exponent 

p = l + e for ^ < 0.5 

p = l-e + e/ii for 0.5 < /!< 2 



(63) 
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for small distances r < D i^/^. This prediction is valid only in the case fj, > 6 for which the series (I61> is 
convergent. However, the same asymptotic results are also obtained by different methods in the case ji < 6, for 
instance expression ( I63> is recovered for all /i < 2 using the integral representation of |89] [A. Saichev, private 
communication]. The numerical evaluation of (I56> . which converges for ji < 6, also recovers the asymptotic 
results (I62> . The two regimes < 0.5 and 0.5 < /i < 2 are illustrated in Figures and |8] respectively. The 
seismicity rate N{t, r) is evaluated from expression (I56> for small z and from expression (I61> for large z. 

We also performed numerical simulations of the ETAS and CTRW models and the results are in good agreement 
with expression j56> and ibH for N{'r, t) for t c and r d. For very small times t <C c, or for very small 
distances r <^ d, expressions i56\ and i6l\ are not valid because they are based on a long wavelength |fc| 
and long time /3 — > approximation. Numerical simulations of the ETAS model in the case 6 — 0.2 and fi — 0.9 
are presented in Figure|9l and are in good agreement with the analytical solutions (15 6> and i61\ shown in Figure|8] 
for the same parameters, except from the truncation of N{t, r) for times t ^ c and distances r <^ d that are not 
reproduced by the analytical solution. 



F. A simple non-separable joint distribution of waiting times and jump sizes: coupled spatial diffusion and long 

waiting time distribution 

Consider the choice for (^,„. [t — ti,r — fi) replacing (|3 by 

<j}mAt-U,f-n)^ p{mi) ^{t~u) <^{\f~n\l^/Dt) , (64) 
where p{mi) and '^{t) are again given by Q and while (jS) is changed into 

1 



$(|r - n\HDt) = ^= exp (-|f - nY/Dt) . (65) 

The spatial diffusion of seismic activity is now coupled to the waiting time distribution. Expression \65\ captures 
the effect that, in order for aftershocks to spread over large distances by the underlying physical process, they need 
time. In fact, returning to the discussion in the introduction on the various proposed mechanisms for aftershocks, 
expression ( I65> embodies a microscopic diffusion process. 
In this case, \2\\ must be replaced by 



where 0(/3, k) is the Laplace-Fourier transform of the product ^{t) $(|r|/V-Di). For large times and long dis- 
tances for which the first terms in the expansion in (3 and k are sufficient, and for n = 1, we obtain 

The inverse Laplace-Fourier transform of (I66> is 

^^''"^ ^ ^ V^t °^P(-|^1'/^i) ■ (68) 

As expected, expression (I68> recovers the dressed Omori propagator in the case of absence of space dependence 
fid]. At finite r and long times, the dressed Omori law also decay as l/t^~^. The diffusion of aftershocks is 
normal with the standard diffusion exponent H — 1/2. 



V. NEW QUESTIONS ON AFTERSHOCKS DERIVED FROM THE CTRW ANALOGY 

We list a series of comments and questions suggested from the analogy between the ETAS model and the CTRW 
model. In particular, we discuss the possibility of defining new observables for earthquake aftershocks, that could 
be worthwhile to investigate in future empirical studies of earthquake aftershocks. 
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A. Recurrence of aftershock activity in tlie proximity of the main shock 

A quantity often investigated in studies of random walks is the probability W{t, 0) to find the random walker at 
its starting point (the origin) at time t. In the earthquake framework, this is the seismic aftershock rate close to the 
main shock. 



B. First-passage times 

The first passage time of a random walk is the first arrival time of the random walk at a given point r. In 
the earthquake context, this translates into the study of the waiting time for a given region to have its own first 
aftershock after the main shock occurs. The distribution of such first passage waiting times gives the distribution 
of times with no nearby seismic activity. See for instance 185] in the case of a power law distribution of waiting 
times and Gaussian distribution of jump sizes. Margolin and Berkowitz |77] give the distribution of first-passage 
times in the case where the jump distribution is narrow and the waiting distribution is long-tailed ~ They 
analyze the three different regimes 6 < 1, 1 < 6 < 2 and 6 > 2. 

C. Occupation time of seismic activity 

Weiss and Calabrese lloill have studied the total amount of time spent by a lattice CTRW on a subset of points. 
In the seismic language, this amounts to study the probability distribution of the durations of aftershock sequences 
that are localized in a specific subset of the space. In other words, how probable are aftershock sequences that are 
found only within a given spatial subset over a certain duration? 

D. Transience and recurrence of seismic activity 

Another question that has been studied in some details in the CTRW framework is whether random walks are 
transient or recurrent. A transient random walk visits any point r at most a finite number of times before escaping to 
infinity. For earthquakes, the transient regime corresponds to the activation of at most a finite number of aftershocks 
in any given point r. In contrast, a recurrent random walk may return a growing number of times to all or a subset 
of points at time increases. In the aftershock language, this means that these points wiU have a never-ending 
(decaying) aftershock activity. We stress here the difference between the global Omori law giving a never-ending 
power law decay of the aftershock activity (in the sub-critical regime n < 1) and its spatial dependence which must 
exhibit important variations. In particular, in the recurrent regime, an Omori law can be documented by counting 
aftershocks in those limited regions of space which are activated again and again. 

E. Probability for the cumulative number of aftershocks 

Let us define a basic quantity in the CTRW formalism, namely the probability Xrn{t) to make exactly m steps up 
to time t. In the earthquake context, Xm (t) is the probability to have exactly m aftershocks after the main shock. In 
the case in which the spatial transition probability i'(r) between different positions is independent of the waiting 
times (corresponding to factorizing {t — ti, f— -fi) as in (|3), the probability density W{t, r) to find the walker 
at position fat time t can be written 

W{t,r) = Y,Wm{T^Xm{t) , (69) 

m=0 

where Wmir) is the probability to reach f from in m steps. In the earthquake context, Wmir) is the probability 
that there has been exactly m events in the time interval [0, t] and that the last one occurred at r. Equation ( I69> 
states that the CTRW is a random process subordinated to simple random walks described by Wm (f) under the 
operational time given by the Xm{t) distribution I92ll93ll . 
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F. Random walk models with birth and death and background seismicity from localized sources 

Bender et al. have studied models of random walks in which walkers are born in proportion to the population 
at one specific site (for instance the origin) with probability a — 1 (with a > 1) and die at all other sites with 
probability 1 — n (with n < 1). In the earthquake context, this consists in assuming that the aftershock activity is 
fed by a localized region in space, which is itself activated by the aftershocks returning to this region, furthering 
the overall activity. This may be considered to describe the seismic activity close to a plate boundary, in which 
the plate boundary is the constant self-consistent source of a seismic activity which may spread over a significant 
region away from the boundary. The excursion of the random walkers quantify the spread of the seismic activity 
away from the main fault structure. The rate of death of the walkers correspond exactly to the distance 1 — n from 
the critical value n = 1. Bender et al. |94] find a phase diagram in the (a — 1, 1 — n) parameter space in which a 
boundary separates two possible asymptotic regimes: 

1 . for small a — 1 and large 1 — n, the seismic activity at the origin and everywhere eventually dies off; 

2. for large a — 1 and small 1 — n, the average seismic activity at the origin approaches a positive constant at 
long times. In this regime, there is a transition as a — 1 is decreased or as 1 — n is increased, between a case 
where the global seismic activity outside the origin goes to zero and a case where it diverges at long times. 
On the boundary between these two regimes in the (a — 1 , 1 — n) parameter space, the distribution of seismic 
activity approaches a steady state at long times. There is a critical point (for space dimensions different from 
2) at a certain value (oc — 1, 1 — ric), for which the long-time seismic activity away from the source is given 
by ~ (a — flc)" where is a critical exponent equal to 2 in three dimensions. 

Note that the results of |94] are obtained for random walks on a lattice. This can easily be converted into a CTRW 
by the fact that a CTRW is nothing by a process subordinated to discrete random walks under the operational time 
defined by the process {ti} of the time of just arrival to a given site, as given by (I69> . 

VI. DISCUSSION 

Using the analogy between the ETAS model and the CTRW model established here, we have derived the relation 
between the average distance between aftershocks and the mainshock as a function of the time from the mainshock, 
and the joint probability distribution of the times and locations of aftershocks. 

We have assumed that each earthquake triggers aftershocks at a distance r and time t according to the bare 
propagator (j){r, t), which can be factorized as 5'(t)$(r). This means that the distribution $(r) of the distances be- 
tween an event and its direct aftershocks is decoupled from the distribution {t) of waiting time. Hence, the direct 
aftershocks triggered by a single mainshock do not diffuse in space with time. Notwithstanding this decoupling 
in space and time of the bare propagator (f){r, t), we have shown that the global law or dressed propagator N{t, r) 
defined as the global rate of events at time t and at position r, cannot be factorized into two distributions of waiting 
times and space jumps. This joint distribution of waiting times and positions of the whole sequence of aftershocks 
cascading from a mainshock is different from the product of the bare time and space propagators. 

The mean distance between the mainshock and its aftershocks, including secondary aftershocks, increases with 
the time from the mainshock, due to the cascade process of aftershocks triggering aftershocks triggering after- 
shocks, and so on. In the critical case n = 1, this diffusion takes the form of a power-law relation R ^ of 
the average distance R between aftershocks and the mainshock, as a function of the time t from the mainshock. 
If the local Omori law is characterized by an exponent < 6* < 1, and if the space jumps follow a power law 
$(r) ^ l/(r + (i)^+^, the diffusion exponent is given by iJ = 6/ ji'm the case n < 2 and H — 9/2 in the case 
ji > 2. Depending on the 9 and ji values, we can thus observe either sub-diffusion (H < 1/2) or super-diffusion 
(H > 1/2), as summarized in Figure [TO] In the sub-critical (n < 1) and super-critical (n > 1) regimes, this 
relation is still valid up to the characteristic time t* given by Q and for distances smaller than r* c>c Dt*^ given 
by ( 13 0> . For t > t* and r > r* in the sub-critical regime, the global distributions of times and distances between 
the mainshock and its aftershocks are decoupled and there is therefore no diffusion. In the super-critical regime, 
the aftershock rate increases exponentially for t > t* and the aftershocks diffuses more rapidly than before t*. 

In the critical regime, the cascade of secondary aftershocks introduces a variation of the apparent Omori ex- 
ponent as a function of the distance from the mainshock. The asymptotic values of the Omori exponent in the 
different regimes are summarized in Table HII In the regime /i < 2, we observe a transition from an Omori law 
decay with an exponent p = 1 — 26* at early times ^ r/D to a larger exponent at large times. This provides 
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another mechanism to explain the observed variability of the Omori exponent. In the regime /i > 2, a power-law 
decay of the seismicity with time is observed only at large times ^ r/D. At early times, or at large distances 
r 3> Dt^ , the seismicity rate is very small, because the seismicity as not yet diffused up to the distance r. 

We should emphasize that our theoretical analysis of aftershock diffusion predicts the behavior of the ensemble 
average of aftershock sequences. Individual sequences may depart from this ensemble average, especially for 
sequences with few earthquakes and limited durations. For long sequences (20,000 events say), we have verified 
that the exponent H measured on individual sequences does not deviate from the ensemble average value by more 
than about 20%. As already discussed, the impact of fluctuations become however more effective as the parameter 
a increases above 6/2. 

The diffusion of the seismicity also renormalizes the spatial distribution of the seismicity, which is very different 
from the local distribution ^{r) of distances between a triggering event and its direct aftershocks. In the regime 
/i > 2, the global seismicity rate N{t, f) decays exponentially with the distance from the mainshock, whereas the 
local distribution of distances <l>(r) is a power-law distribution. In the regime /i < 2, the local law $(r) ^ r^^^^ 
is recovered at large distances, but a slower decay for 0.5 < /i < 2 or a constant rate for fi < 0.5 is observed 
at small distances r ^ Dt^ . These predictions on the decrease of the Omori exponent with r have not yet 
been observed in earthquake catalogs, but an expansion of the aftershock zone has been reported in many studies 
128. .22, ,30^ 2Zi 22i 2^ 25^ .36.1 . However, very few studies have quantified the diffusion law. Noir et al. 
I3j1 show that the earthquake Dobi sequence (central Afar, August 1989) composed of 22 M > 4.6 earthquakes 
presented a migration that was in agreement with a diffusion process due to fluid transfer in the crust, characterized 
by a normal diffusion process with exponent H = 0.5. Tajima and Kanamori |31, 32] studied several aftershock 
sequences in subduction zone and observed a much slower logarithmic diffusion, which is compatible with a low 
exponent H close to 0.1. In some cases, the aftershock sequence displays no expansion with time. For instance, 
Shaw 1 95] studied several aftershock sequences in California and concluded that the distribution of distances 
between the mainshock and its aftershocks is independent of time. This can be explained by the fact that the Omori 
exponent measured in L95.1 is very close to 1, thus 9 is very small and our prediction is that the exponent H should 
be very small. 

In fact, the ETAS model predicts that diffusion should be observed only for aftershock sequences with a mea- 
sured Omori exponent p significantly smaller than 1, which can only occur according to our model when the bare 
Omori propagator with exponent 1 + 6 is renormalized into the dressed propagator with global exponent 1 — 9. We 
have shown that this renormalization of the exponent only occurs at times less than t*, while for longer times in 
the sub-critical regime n < 1 the dressed Omori propagator recovers the value of the bare exponent 1 + 6 > 1 (see 
figure|2}. Therefore, identifying an empirical observation of p < 1 with our prediction p =1 — 6 indicates that the 
aftershock sequence falls in the "good" time window t < t* in which the renormalization operates. We have also 
shown that the dressed propagator gives a diffusion only for t < t*. We can thus conclude that, according to the 
ETAS model, the observation of an empirical Omori exponent larger than 1 is indicative of the large time t > t* 
behavior in the sub-critical regime n < 1^ for which there is no diffusion. This provides a possible explanation for 
why many sequences studied by IBU l32L 19511 do not show a diffusion of the aftershock epicenters. Reciprocally, a 
prerequisite for observing diffusion in a given aftershock sequence is that the empirical p-value be less than 1 in 
order to qualify the regime t < t*. 

An alternative model has been discussed by Dieterich ll40ll who showed that the spatial variability of the stress 
induced by a mainshock, coupled with a rate and state friction law, results in an expansion of the aftershock zone 
with time. This expansion does not take the form of a diffusion law as observed in the ETAS model, the relation 
between the characteristic size of the aftershock zone does not grow as a power law of the time from the mainshock 
(equation (22) and Figure 6 of ll40ll '). 

Marsan et al. ]96, 97] and Marsan and Bean ig^l studied several catalogs at different scales, from the scale of a 
deep mine to the world-wide seismicity, and observed that the average distance between two earthquakes increases 
as a power-law of the time between them, with an exponent often close to 0.2, indicative of a sub-diffusion process. 
They interpreted their results as a mechanism of stress diffusion, that may be due to fluid transfer with heteroge- 
neous permeability leading to sub-diffusion. Their analysis is quite different from those used in other studies, 
because they consider all pairs of events, without distinction between aftershocks and mainshocks. This analysis 
can however lead to spurious diffusion, and in some cases this method does not detect diffusion in synthetic data 
set with genuine diffusion. We have tested their analysis on a synthetic catalog generated by superposing a back- 
ground seismicity with uniform spatial and temporal distribution, and 10 mainshocks with poissonian distribution 
in time and space, and with a power-law distribution of energies. Each of these mainshocks generates only direct 
aftershocks, without secondary cascades of aftershocks, and the number of aftershocks increases exponentially 
with the magnitude of the mainshock. This way, we generate a synthetic catalog without any physical process of 
diffusion, and which includes all the other well-established characteristics of real seismicity: clustering in space 



23 



and time superposed to a seismicity background. Applying the analysis of to this synthetic data set 

leads to an apparent diffusion process with a well-defined exponent H — 0.5. However, this apparent diffusion 
does not reflect a genuine diffusion but simply describes the crossover from the characteristic size of an aftershock 
zone at early times to the larger average distance between uncorrected events at large times. In plain words, the 
apparent power law R oc is nothing but a cross-over and is not real. Furthermore, applying this analysis to a 
synthetic catalog generated using the ETAS model, without seismicity background, and with a theoretical diffusion 
exponent H = 0.2, the method yields H = 0.01 if we use all the events of the catalog. If we select only events up 
to a maximum distance r„iax to apply the same procedure as in f9^ l97ll98il . we obtain larger values of H which 
are more in agreement with the theoretical exponent H = 0.2 but with large fluctuations that are function of r^ax- 
Therefore, it is probable that the diffusion reported in |96, is not real and results from a cross-over between 

two characteristic scales of the spatial earthquake distribution. It may be attributed to the analyzing methodology 
which mixes up uncorrelated events. We are thus reluctant to compare the results of Marsan et al. |96,97,98] with 
the predictions obtained with the ETAS model. 

One can similarly question the results on anomalous diffusion of seismicity obtained by Sotolongo-Costa et 
al. i99il . who considered 7500 micro-earthquakes recorded by a local Spanish network from 1985 to 1995. They 
interpret the sequence of earthquakes as a random walk process, in which the walker jumps from an earthquake 
epicenter to the next in sequential order The time between two successive events is seen as a waiting time between 
two jumps and the distance between these events is taken to correspond to the jump size. Since the distributions of 
time intervals and of distances between successive earthquakes are both heavy-tailed (approximately power laws), 
their model is a CTRW. We cannot stress enough that their CTRW model of seismicity has nothing to do with 
our results on the mapping of the ETAS model onto a CTRW. Their procedure is ad-hoc and their results depend 
obviously strongly on the space domain of the analysis since distant earthquakes that are completely unrelated can 
be almost simultaneous! We also stress that our mapping of the ETAS model onto the CTRW model does not 
correspond to identifying_^an earthquake sequence as a single realization of a CTRW, as assumed arbitrarily by 
Sotolongo-Costa et al. \9%. 

Our predictions obtained here are thus difficult to test on seismicity data, due to the small number of events 
available and the restricted time periods and distance ranges, and because the seismicity background can strongly 
affect the results. New methods should hence be developed to investigate if there is a real physical process of 
diffusion in seismic activity and to compare the observations of real seismicity with the quantitative predictions 
of the ETAS model. Preliminary study of aftershock sequences in California leads to the conclusion that most 
aftershock sequences are characterized by an Omori exponent p > 1, indicative of the sub-critical regime with 
t > t* . As expected from our predictions in this regime, we do not observe an expansion of the aftershock zone. 
However, a few sequences give a value p < 1 and also exhibit an increase of the average distance between the 
mainshock and its aftershocks consistent with our predictions. A detailed report of this analysis will be reported 
elsewhere. 



VII. CONCLUSION 

We have studied analytically and numerically the ETAS (epidemic-type aftershock) model, which is a simple 
stochastic process modeling seismicity, based on the two best-established empirical laws for earthquakes, the 
power law decay of seismicity after an earthquake and a power law distribution of earthquake energies. This 
model assumes that each earthquake can trigger aftershocks, with a rate increasing with its magnitude. In this 
model, the seismicity rate is the result of the whole cascade of direct and secondary aftershocks. 

We have first established an exact correspondence between the ETAS model and the CTRW (continuous-time 
random walk) model. We have then used this analogy to derive the joint probability of times and distances of the 
seismicity following a large earthquake and we have characterized the different regimes of diffusion. 

We have shown that the diffusion of the seismicity should be observed only for times t < t*, where t* is a 
characteristic time depending on the model parameters, corresponding to an observed Omori exponent smaller 
than one. Most aftershock sequences have an observed Omori exponent larger than one, corresponding to the 
subcritical regime of the ETAS model, for which there is no diffusion. The diffusion of the seismicity produces a 
decrease of the Omori exponent as a function of the distance from the mainshock, the decay of aftershocks being 
faster close to the mainshock than at large distances. The spatial distribution of seismicity is also renormalized by 
the cascade process, so that the observed distribution of distances between the mainshock and its aftershocks can 
be fundamentally different from the bare propagator $(r) which gives the distribution of the distances between 
triggered and triggering earthquakes. We have also noted that the ETAS model generates apparent but realistical 
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fractal spatial patterns. 

Assuming that the distances between triggering and triggered events are independent of the time between them, 
this model generates a diffusion of the whole sequence of aftershocks with the time from the mainshock, which 
is induced by the cascade of aftershocks triggering aftershocks, and so on. Our results thus provides a simple 
explanation of the diffusion of aftershock sequences reported by several studies, which was often interpreted as 
a mechanism of anomalous stress diffusion. We see that no such "anomalous stress diffusion" is needed and our 
theory provides a parsimonious account of aftershock diffusion resulting from the minimum physical ingredients 
of the ETAS model. As Einstein once said, "A theory is more impressive the greater the simplicity of its premises, 
the more different the kinds of things it relates and the more extended its range of applicabihty." 
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ETAS 


CTRW 




pdf for a "daughter" to be born at time t 
from the mother that was born at time 


pdf of waiting times 


$(r) 


pdf for a daughter to be triggered 
at a distance r from its mother 


pdf of jump sizes 


m 


earthquake magnitude 


tag associated with each jump 


p(m) 


number of daughters 
per mother of magnitude m 


local branching ratio 


n 


average number of daughters created per mother 
summed over all possible magnitudes 


control parameter of the random 
walk survival (branching ratio) 


n <1 


subcritical aftershock regime 


subcritical "birth and death" 


n = 1 


critical aftershock regime 


the standard CTRW 


71 > 1 


supercritical exponentially 
growing regime 


explosive regime of the 
"birth and death" CTRW 




number of events of any possible 
magnitude at r at time t 


pdf of just having 
arrived at f at time t 


W{t,f) 


pdf that an event at r has occurred at a time t' < t 
and that no event occurred anywhere from t' to £ 


pdf of being at r at time t 



TABLE I: Correspondence between the ETAS (Epidemic-type aftershock sequence) and CTRW (continuous-time random 
walk) models, 'pdf stands for probability density function. 





large z 


small z 




r ^Dt" 




H < 0.5 


p = i + e 


p = l-2e 


0.5 < < 2 


p = i-e + e/fi 


p = l-2e 


2 < M 


p =1-6/2 


not defined " 



"The Omori exponent is not defined in this case because the dependence of N(t,r) with respect to time given by expression l42l and 
represented in Figurel4lhas a contribution from the exponential asymptotics which is different from a power-law for large distances r ^ D . 

TABLE IL Asymptotic values of the (renormalized) Omori exponent (of the dressed propagator) in the different regimes for 
z <C 1 and z S> 1 where z = S-L- . 
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FIG. 1: Maps of seismicity generated by the ETAS model with parameters b = 1, 6 — 0.2, jj, — 1, d = 1 km, a = 0.5, 
c = 0.001 day and a branching ratio n = 1. The mainshock occurs at the origin of space with magnitude Af — 7. The 
minimum magnitude is fixed at mo = 0. The distances between mainshock and aftershocks follow a power-law with parameter 
/i = 1 and the local (or bare) Omori law is oc l/f^"*"". According to the theory developed in the text, the average distance 
between the first mainshock and the aftershocks is thus expected to grow as i? ~ with H = 0.2 (equation i53\ ). The two 
plots are for different time periods of the same numerical simulation, such that the same number of earthquakes A'^ = 3000 
is obtained for each graph: (a) time between and 0.3 days; (b) time between 30 and 70 yrs. Real aftershock sequences are 
indeed observed to last decades up to a century. Large black dots indicate large aftershocks around which other secondary 
aftershocks cluster. The mainshock is shown by a black star. At early times, aftershocks are localized close to the mainshock, 
and then diffuse and cluster around the largest aftershocks. 
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FIG. 2: Seismicity rate N{t) for the temporal ETAS model calculated for 6 — 0.3 and c = 0.001 day. The local law 
(p{t) oc which gives the probability distribution of times between an event and its (first-generation) aftershocks is 

shown as a dashed line. The global law N{t), which includes all secondary and successive aftershocks generated by all the 
aftershocks of the first event, is shown as a solid line for the three regimes, n < 1, n = 1 and n > 1. In the critical regime 
n = 1, the seismicity rate follows a renormalized or dressed Omori law oc l/i*" for t > c with an exponent p = 1 — 0, smaller 
than the exponent of the local law 1 + 6. In the sub-critical regime (n < 1), there is a crossover from an Omori law l/t^~^ 
for t < t* to for t > t* . In the super critical regime (n > 1), there is a crossover from an Omori law l/t^~^ for 

t < t* Xo an exponential increase N{t) ^ exp{t/t*) for t > t* . We have chosen on purpose values of n = 0.9997 < 1 and 
n = 1.0003 > 1 very close to 1 such that the crossover time t* = 10® days given by is very large. In real data, such large 
t* would be undistinguishable from an infinite value corresponding to the critical regime n — 1. This representation is chosen 
for pedagogical purpose to make clear the different regimes occurring at times smaller and larger than t* . In reality, we can 
expect n to be significantly smaller or larger than 1, such that t* becomes maybe of the order of months, years to decades and 
the observed Omori law will thus lie in the cross-over regime, given an apparent Omori exponent anywhere from 1 — 6 to 1 + ^. 
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FIG. 3: Plot of the correlation function of the 3.000 epicenters generated in the time interval [30, 70] yrs and shown in the right 
panel of figureQ calculated following Grassberger-Procaccia's algorithm j^, as a function of scale r, in double-logarithmic 
scales. 
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FIG. 4: Rate of seismicity N(t, r) in the critical regime n = 1 foi 9 — 0.2, /j. > 2, c' — 1 day and a — 1 km, evaluated from 
expressions <40t and <42t . plotted as a function of the time (a) for different values of the distance r between the mainshock 
and its aftershocks, and (b,c) as a function of r (logarithmic scale for r in (b) and linear scale for r in (c)) for different 
values of the time between the mainshock and its aftershocks. The temporal decay of seismicity with time is characterized 
by a power-law decay N{r,t) ~ l/t^^^^'^ close to the mainshock epicenter or at large times for r <^ Dt^^^. For large 
distances r ^ Dt^^^, there is a truncation of the power-law decay at early times t^^^ <^ r/D, because the seismicity has 
not yet diffused up to the distance r. Although the distribution of distances between a mainshock and its direct aftershocks 
$(r) follows a power-law distribution with exponent 1 + /i, the log-linear graph (c) shows that the global rate of aftershocks 
N{f, t) decreases approximately exponentially as a function of the distance from the mainshock, with a characteristic distance 

r 2 

that increases with time. This is in agreement with expression <42t which predicts N{t, r) ~ exp (^|r|/Di^/^) ^"^ 

N(t, r) ~ exp (C(t)lrl'') with an exponent q = 2/(2 — 6) close to 1 within the exponential. The same remark as for figure 
|2| applies: the representation of our predictions for very large times is made for pedagogical purpose to illustrate clearly the 
different regimes. 
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FIG. 5: Average distance between the first mainshoclc and its aftersiiocks as a function of the time from the mainshock, for 
numerical simulations of the ETAS model in the critical regime n = 1, generated with the parameters 9 = 0.2, d = 1 km, 
fj, = 3 and c = 10^'^ day. The theoretical prediction for the diffusion exponent is thus H — 9/2 — 0.1. We observe a crossover 
from a larger exponent at early times when the mean distance is close to the characteristic scale d = 1 km of the distribution 
of distances between an aftershock and its progenitor, to a sub-diffusion with an exponent close to the theoretical prediction at 
large times. The solid line is a fit of the numerical data for times t > 10 days, which gives an exponent H = 0.12 slightly 
larger than the predicted value H — 0.1. 
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FIG. 6: Average distance between the first mainshock and its aftershocks as a function of the time from the mainshock, for a 
numerical simulation of the ETAS model in the critical regime n — 1, with 9 — 0.2, fi = 0.9, c' = 1 day and d = 1 km. The 
solid line is a fit of the data which gives an exponent H — 0.25 in good agreement with the predicted value H = 0.22. 
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FIG. 7: Rate of seismicity N{t, r) for 9 — 0.2, /j, = 0.2, c' — 1 day and a — 1 km, evaluated from expressions J56> and <62t . 
plotted as a function of the time (a) for different values of the distance r between the mainshock and its aftershocks, and (b) as 
a function of r for different values of the time between the mainshock and its aftershocks. We stress again that the time scales 
shown here do not necessarily correspond to real observable time scales but are presented to demonstrate clearly the existence 
of the two regimes. The dashed lines give the predicted asymptotic dependence in each regime. 
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FIG. 8: Rate of seismicity N(t, r) for Q — 0.2, [i = 0.9, d — \ day and a — \ km, evaluated from expressions j56> and <62t . 
plotted as a function of the time (a) for diffe rent values of the distance r between the mainshock and its aftershocks, and (b) as 
a function of r for different values of the time between the mainshock and its aftershocks. The dashed lines give the predicted 
asymptotic dependence in each regime. 
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FIG. 9: Rate of seismicity N{t, r) obtained from numerical simulations of the ETAS model generated with the same parameters 
as in Figure|8|(S — 0.2, fi — 0.9, c' — 1 day and d = 1 km). N(r, t) is computed by averaging over 500 numerical realizations 
of the ETAS model, (a) aftershock rate as a function of the time from the mainshock for several distances |r| ranging from 
0.01 to 10* km. (b) Apparent Omori exponent measured for times t > 10 as a function of the distance from the mainshock. 
The aftershock decay rate (with time) is larger close to the mainshock epicenter than at large distances from the mainshock. 
The asymptotic values for small and large distances are in agreement with the predictions <63> for r <C Dt^^^ and <59t for 
r > Dt"^", which are shown as the horizontal dashed lines, (c) Rate of seismicity N{t, r) as a function of the distance 
between aftershocks and mainshock for various times. The theoretical prediction for large distances is shown as the dashed line 
with slope —(1 + jj,). 
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FIG. 10: Classification of the different regime of the diffusion of aftershocks in space as a function of time from the main 

shock. The bare Omori law for aftershocks decay with time as 1/t^^*. The jump size distribution between the earthquake 
"mother" and its "daughters" is proportional to 1/r^"'"''. R{t) is the average distance between all aftershocks triggered up to 
time t after the mainshock. 



